Critical properties of Ising model on Sierpinski fractals. 
A finite size scaling analysis approach. 



Jose M. Carmona a , Umberto Marini Bettolo Marconi 5 , 
Juan J. Ruiz-Lorenzo c and Alfonso Tarancon" 

a Departamento de Fisica Teorica, Universidad de Zaragoza, 
Plaza S. Francisco s/n, 50009 Zaragoza (Spain) 

5 Dipartimento di Matematica e Fisica and Istituto Nazionale 
di Fisica della Materia, Universita di Camerino, 
Via Madonna delle Carceri, 62032, Camerino (Italy) 

c Dipartimento di Fisica and Istituto Nazionale 
di Fisica Nucleare (Sez. Roma-I), Universita di Roma La Sapienza, 
P. A. Moro 2, 00185 Roma (Italy) 

carmona , tarancon@sol . unizar . es 
umberto . marini . bettoloOromal . inf n . it 
ruizOchimera . romal . inf n . it 

February 1, 2008 

Abstract 

The present paper focuses on the order-disorder transition of an Ising model on 
a self-similar lattice. We present a detailed numerical study, based on the Monte 
Carlo method in conjunction with the finite size scaling method, of the critical prop- 
erties of the Ising model on some two dimensional deterministic fractal lattices with 
different Hausdorff dimensions. Those with finite ramification order do not display 
ordered phases at any finite temperature, whereas the lattices with infinite connectiv- 
ity show genuine critical behavior. In particular we considered two Sierpinski carpets 
constructed using different generators and characterized by Hausdorff dimensions 
d H = In 8/ In 3 = 1.8927.. and d H = In 12/ In 4 = 1.7924.., respectively. The data 
show in a clear way the existence of an order-disorder transition at finite tempera- 
ture in both Sierpinski carpets. By performing several Monte Carlo simulations at 
different temperatures and on lattices of increasing size in conjunction with a finite 
size scaling analysis, we were able to determine numerically the critical exponents 
in each case and to provide an estimate of their errors. Finally we considered the 
hyperscaling relation and found indications that it holds, if one assumes that the 
relevant dimension in this case is the Hausdorff dimension of the lattice. 



1 Introduction 



The present understanding of phase transitions has greatly benefited from the study of 
the spin-lattice models, perhaps the simplest examples of extended systems showing non 
trivial cooperative behavior, such as spontaneously symmetry breaking. In most cases 
one is interested in studying systems whose geometrical properties are regular so that one 
assumes that the spins occupy the cells of a regular Bravais lattices. Near the critical point, 
i.e. when the correlation length is much larger than the lattice spacing, the influence of 
the lattice structure becomes negligible and only the embedding dimension, the number 
of components of the order parameter together with its symmetry, and the nature of the 
couplings concur to determine the values of the critical exponents. 

Such universal behavior is absent if the lattice is a fractal, because the translational 
invariance is replaced by the much weaker dilation invariance. Notwithstanding the intense 
activity on various physical problems in a space which instead of being Euclidean is a 
fractal lattice, the issue of the phase transitions and of the critical properties on self- 
similar supports has been rarely addressed ||, [| ■ Unlike critical phenomena in spaces of 
integer dimension those occurring in self-similar geometries have not been explored so far, 
apart from some isolated cases. One expects that the lack of translational invariance plays 
a crucial role in phase transitions on fractal supports. These systems besides serving in 
practice to model natural materials such as porous rocks, aerogels, sponges etc., provide a 
geometric realization of non integer Hausdorff dimension. Therefore they offer a possibility 
of testing the e-expansion technique for e not integer. Finally, one can explore systems 
whose geometrical dimension is very near to its lower critical dimension (in the Ising Model 
is one: i.e. one is the larger (integer) dimension in which the system does not have a phase 
transition at finite temperature). 

In the present study we shall consider lattices obtained by removing sites from a square 
lattice. If the diluted lattice is generated by a sequence of random deletions one obtains 
the so-called site diluted Ising model (SDIM), whose phase diagram has been studied (see 
for instance @, [5|, ||). In this case, one finds that the critical temperature T c tends to zero 
as the probability p of having a site tends to the percolation threshold p c = 0.592746 (on a 
square lattice). Notice that only at the percolation threshold p c the lattice manifests true 
self-similarity (and fractal dimension 1.8958 0), but the phase transition occurs only at 
T = 0. Thus to observe simultaneously genuine criticality and self-similarity we consider 
a non-stochastic lattice, the so called Sierpinski carpet (SC). As pointed out some years 
ago, one can decide whether a lattice is able to support an order-disorder transition of the 
Ising type by looking at its order of ramification TZ, which is finite if after eliminating a 
finite number of bonds one can isolate an arbitrarily large sublattice. Only if TZ = oo a 
phase transition occurs. Thus the SC displays a critical point at finite temperature, while 
the Sierpinski gasket does not, being characterized by a finite TZ. This is a manifestation 
of the fact that, when considering non regular structures, the embedding dimension alone 
is not sufficient to determine the phase behavior. 

The deterministic SC we consider, in spite of having a lower Hausdorff dimension than 
that of the SDIM, is always above the percolation threshold by construction. One can thus 
observe the interplay between fractal geometry and thermodynamics. 

In the present work we shall investigate the critical behavior of the Ising model on 
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two classes of fractal lattices, the two dimensional Sierpinski gasket and the two dimen- 
sional Sierpinski carpet. The first lattice, with Hausdorff dimension du = In 15/ In 5 = 
1.682606..., can be studied analytically and an exact real space renormalization group 
treatment rules out the possibility of a finite ordering temperature. The Sierpinski car- 
pet instead displays a genuine transition at finite temperature that we have characterized 
numerically for two different fractal dimensions, du = In8/ln3, and du = Inl2/ln4. 

In statistical mechanics it is usually assumed that in the infinite volume limit (V — > oo) 
the average of a given observable calculated over a subvolume sufficiently large compared 
with the bulk correlation length yields a result independent of the particular choice and 
location of the cell within the sample and moreover that this average converges to the 
average value over the whole sample. In the fractal lattices that we consider such property 
breaks down below the critical point because the system cannot be regarded as a uniform 
material over length scales larger than the correlation length, due to the presence of voids 
of all sizes which render the system effectively not uniform. Even the correlation length 
depends on the position. However we shall show that this system has a second order phase 
transition and that means the existence of a divergent length scale ^ in the thermodynamic 
limit. Although the fractal is not homogeneous, we postulate the existence of this "average" 
correlation length at least near criticality. To assume the existence of such a length is in 
some sense similar to what happens when studying the diffusion problem on a lattice. 
There, of course, a walker does not diffuse with the same law from any point of the lattice, 
but one can still define a diffusive type of behavior of the type R ~ t 1 l dw by averaging over 
all possible initial positions of the walker. 

We outline the plan of the paper: in section 2 we introduce the lattices and define 
the observables of the Ising model which will be measured in the simulations, and give 
some details about the numerical procedure employed to obtain the statistical averages. In 
section 3 we illustrate some analytical results concerning lattices with finite ramification 
order, which were obtained by means of an exact Migdal-Kadanoff decimation procedure. 
In section 4 we recall briefly the statement of the finite size scaling method. In section 5 we 
illustrate the results of the simulations for the fractal characterized by du = In 8/ In 3 and 
discuss at some length the details of the data analysis on the basis of which we determined 
the values of the various critical exponents and of the critical coupling. In section 6 we 
illustrate the results of the simulations for the fractal with du = In 12/ In 4, with special 
emphasis to the hyperscaling relation. Finally in 7 we present the conclusions. 

2 The Model and Observables 

Let us introduce the lattice models that we shall consider in numerical studies. We have 
first considered a SC, named hereafter fractal A, constructed starting from a square lattice 
of L x L cells with L = 3 n , dividing it into 3x3 blocks of equal size and discarding the 
cells contained in the central block. Divide again each of the remaining blocks into 3x3 
sub-blocks and discard all the central elements. Carry on this procedure until the smallest 
sub-block contains a single cell. The resulting structure is formed by V cells, where V is 
related to the linear dimension through the Hausdorff dimension du as V = L dH . In this 
case, dn = In 8/ In 3 = 1.892789... Of course one can build many different SC with the 



3 



o 



o 



o 



o o o 



o 



o o o 



o 



o o o 



o 



o 



o 



Figure 1: SC of dimension da — In 8/ In 3 for L = 9. Filled circles represent cells to which 
Ising spins are assigned, whereas empty circles represent cells which have been eliminated 
from the original square lattice. The links represent the interactions among spins. 



same fractal dimension, but having different distributions of voids. This can be quantified 
by means of the so called lacunarity ||. In the present work we stick to the symmetrical 
fractal shown in Fig. [I]. We define bonds between the remaining cells. The number of 
bonds N is 

n = \v + \l, (1) 

so that N is proportional to V in the thermodynamic limit. Note that this fractal has 
1Z = oo. 

For the fractal B, the construction follows the same procedure, but now one considers 
4x4 blocks and discards the four central ones. This lattice has fractal dimension dn = 
In 12/ In 3 = 1.792481... and in this case the relation between the number of bonds N, the 
linear dimension L and the number of cells V is 

N=-V + -L. (2) 
2 2 v ; 

Now, at each unit cell we assign a dichotomic spin variable a = ±1. The total energy 
of the system with ferromagnetic interactions is defined as 

H = - a i a 3 ' ( 3 ) 

<i,j> 

where < i,j > are nearest neighbor cells. The energy normalized to the number of links is 
a quantity that has a well-defined thermodynamic limit, and the specific heat is: 

C = N ({E 2 } - (E) 2 ) . (4) 
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Since our primary target is to obtain informations about the critical behavior of the 
model we first locate the critical point and then we extract the critical exponents. To 
achieve this goal we introduce the observables which are relevant. These are defined as 
follows: the intensive magnetization M 

i 

the isothermal susceptibility: 

X = V{(M 2 )-(\M\f) , (6) 

and the Binder cumulant 

/ = 2 V" (M 2 ) 2 

3 Migdal-Kadanoff Method 

The Migdal-Kadanoff (MK) decimation is an approximation very suitable for analyzing 
low dimensional systems. MK is exact in one dimension, but lacks predictability in higher 
dimensions. The starting point is the formula for the variation of the inverse temperature 
(P) with the scale factor in the MK approximation in a system of dimensionality D |9|, |T0f , 

Tt = g{(3h (8) 

g (/3) = (D — l)f3 + sinh(/3) cosh(/5) In [tanh(/3)] , (9) 

where t is the logarithm of the scale factor in the block construction. 

The zeroes of g{{3) yield the inverse critical temperatures of the system. Also, it is 
possible to compute the v exponent from 



1 dg(P) 



v dp 



(10) 



Such a formula is very interesting because one can obtain the behavior of the critical 
exponents and the critical temperature near the lower critical dimension of the Ising model; 
for instance when D — > 1 

Pc - —, r , v ~ — - — . (11) 

Hc 2(D - 1) ' D - 1 V ; 

If one assumes that the dimensionality of the system corresponds to the Hausdorff 
dimension (i.e. D = dn = 1.892789...) one obtains p c ~ 0.56 and v ~ 1.12. Numerically 
(solving g(Pc) = and 1/u = g'(P c ) for D = 1.89, that corresponds to fractal A) one 
obtains: 

p c = 0.5120 , v = 1.409 . (12) 
For the fractal B (D = 1.79) the results are 

p c = 0.5906 , v = 1.511 . (13) 
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Figure 2: Sierpinski gasket (SG). The spins live in the sites of the construction. We have 
marked with /x the spins that are decimated. See the text for more details. 



Obviously this approach does not take into account the ramification and lacunarity 
of the model. Nonetheless in this section we have obtained a first guess of the critical 
temperature and of the v exponent for a system with the same fractal dimension that the 
Sierpinski carpet (in the most symmetric version) that we will study numerically in the 
next sections. 



3.1 Exact solution of the Sierpinski gasket 



The decimation can be carried out exactly on the Sierpinski gasket, recovering the same 
lattice after a decimation transformation. 

We express the Boltzmann weight in the usual way 



exp(/3si,s 2 ) = cosh(/5) [1 + s\s 2 tanh(/?)] 



(14) 



where we have denoted by s\ and s 2 two generic spins belonging to the lattice. We denote 
by /i the spins that we decimate and by a the rest of the spins. 
The decimation over the \x spins is the following (see Fig. H) 

y, exp [/3((7i/Ji + <7i/! 2 + /ii/i 2 + H\(J<1 + /il/i 3 + A*2/^3 + At20"3 + + /U3O2)] 

Ml>A t 2,A t 3=±l 

= Fy](l + fc(7i/ii)(l + fecri/i 2 )(l + £;/ii/U 2 )(l + knia 2 )(l + kfii^z) 
x (1 + /c/i 2 /i 3 ) (1 + fc/^as) (1 + k/J, 3 a 3 ) (1 + /c/i3Cr 2 ) 
= 8(1 + A;) 3 (l + A: 2 )(l — 3/c + 5fc 2 — 3A: 3 + fc 4 ) 
k 2 



1 + 



l-3k + 5k 2 - 3k 3 + k 4 



(15) 
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where we have defined k = tanh(/3) and F = cosh 9 (/3). 

The Sierpinski gasket is self-similar under this decimation transformation and we can 
see that eq. fljjj) corresponds to a nearest neighbor Ising interaction with a renormalized 
coupling (p R ): i.e. 



exp \fi R {a x a 2 + a 2 a 3 + o 3 oi)] = F'(l + k R <Tx a 2 ) (1 + fcfl<7i<7 3 )(l + k R a 2 a^) 



F'(l + k R ) 



k R (l + £;#) 

1 H 1 , 3 (0"l0"2 + <7lO-3 + 020-3) 



(16) 



and we have denoted with k R = tanh(/?#) the renormalized k and F' = cosh 3 (/?#). We 
finally arrive at the following exact recursion (by matching the coefficient of (01 o 2 + 020-3 + 
o 3 o-i) of eqs. (0) and ([16])) 

k 2 k R (l + k R ) 



1 - 3fc + 5A; 2 - 3A; 3 + A; 4 1 + jfej. 

The critical points correspond to the solutions of the previous recursion formula. The 
solutions are k = and k = 1, i.e. are (3 = and (3 = 00: i.e. there is no phase transition 
at any finite temperature. From Fig. ^ it is clear that TZ remains finite even when the 
linear size diverges. 



4 Finite Size Scaling Method (FSS) 

For a finite system whose typical linear size is L, which is assumed to be much larger 
than the lattice spacing a, the finite size scaling hypothesis |Tl|, postulates that upon 
approaching the critical point the average value of a given observable, P, depends on the 
size and on the temperature through the following scaling relation: 

Pi(t) - / f AY (is) 



Poo(t) \U(t). 

where T c and t are respectively the critical and the reduced temperature t = \T — T c \/T c of 
the system, and represents the value of P in the infinite volume limit. If Poo(t) behaves 
as t~ p when t — > 0, since the correlation length diverges as ^ ~ t~ u , it is clear that -Pl(^) 
will saturate when the ^ becomes comparable with L. This can be formalized by writing: 

P L (t) = L^g {tL l '») , (19) 

where the scaling function g(x) has the limiting behavior g(x) = const as x — >• and 
g ~ x~ p as x — ► 00. 

Thus at the critical point one can write the following finite size scaling formulae for the 
observables (see for instance 0): 

X (L,t = 0)ocL^[l + O(L-% (20) 
C{L,t = 0) ocL a/u [1 + OiL-")} , (21) 
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where u is the correction-to-scaling exponent (it is just the slope of the field theoretical 
/3-function at the fixed point or equivalently the largest irrelevant scaling exponent in the 
Wilson renormalization group), and t is the reduced temperature of the model. 

For non dimensional quantities, denoted by A, in the proximity of the critical point 
(t < 1) one has the scalings: 

A{L, t) = f A (L l /»t) + L-g A (L l ^t) + 0{L^). (22) 

where Ja and qa are scaling functions. Their derivatives with respect to (3 behave as: 

H A 

— (L,t = 0)a^[l + O(r)]. (23) 

In particular in the present work we have chosen A to be the Binder cumulant U and 
A = ln(|M|). 

The uj exponent can be calculated using the effective critical exponents that are ex- 
tracted from the peaks of the observables at lattice sizes L and sL. Let us take the 
susceptibility, for example. We obtain the effective %(L,sL) as 

l(L,sL) = ±- ln^. (24) 
v Ins x\L) 

Using the scaling formulae and working in the large lattices, we obtain 

-(L,sL) = - + BL~", (25) 
v v 

where 7/V is the infinite volume extrapolation of the ratio 7/V(L, sL) measured on finite 
lattices. We obtain u and 7/1/ by fitting the data to formula (p5|). 

In order to determine the critical temperature, let us call P®(L) the value of (3 where the 
observable O (with positive dimension) displays a maximum. In order to assess its location 
several simulations at different values of f3 are needed in principle, and even doing so its 
precise determination is hard. However, the spectral density method (SDM) renders 
such task easier and much faster. In fact, by using the data from a single simulation at a 
given temperature, say /3q, one can obtain information about some (3\ in the neighbourhood 
of /3q. With SDM one usually gains one order of magnitude in the accuracy of the location 
of the peak. 

The idea is to write the partition function under the form: 

Z oc e ~ m{a) = EE e ~ ma) 6 ( H - E ) = Y. e- 0E N(E), (26) 

{confs} E {confs} E 

where TC(cr) is the Hamiltonian as a function of the configuration {a} and N(E) is the 
energy density of states, which can be obtained from the number of times that the value 
E of the energy is generated during a Monte Carlo simulation. Given the density of states 
at (3 = /3 , the value of a certain operator O at Pi ^ p is given by 

(0(g . _ Y. E O^E)N{E)e-^ E 

{0{(3l)) ~ Z E N(E)e-<*-W • (27) 
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This extrapolation is exact. The problem is that for large deviations A/3 = (3\ — (3q the 
errors in the extrapolation are very large, and then one has to restrict oneself to small A/3 
values (of order 1/y/V). 

In d > 2 the value /3^(L) hardly changes with O, and with the help of the SDM one 
can reach the peak of every observable with a single simulation at a certain (3. We found 
that this is not the case on fractal lattices with du < 2 because the transition region turns 
out to be very large, and sometimes we were forced to perform different simulations at 
different points, depending on the observables, in order to get their corresponding peaks. 

To compute thermal averages we have employed a combination of m steps of the classical 



Metropolis algorithm followed by n steps of the Wolff single-cluster algorithm pA \, which 



gives a very short autocorrelation time. We checked that such method is sufficient to 
ensure ergodicity. However, the Wolff algorithm is unable to flip clusters of intermediate 
size, which might be important in our model, where there are domains of spins at every 



scale. That is why we have applied in addition a Swendsen-Wang algorithm [15|, which is 



more consuming in computing time than the Wolff algorithm, but is able to flip domains of 
spins of every scale. In fact, we have compared the results of the two different procedures. 
Fig. [3] shows that the two methods give compatible results in a region around the simulation 
point, where the SDM extrapolation is valid (away from this region, the SDM extrapolation 
gives large errors and it is therefore not reliable, so that it is not surprising that it gives 
different results for the two algorithms). We also note the narrow range of validity of the 
SDM which is used to extrapolate the numerical data, which makes even more troublesome 
the problem of finding the peaks of observables in order to perform FSS. As far as the 
observables directly measured are concerned one can hardly see differences in the two 
simulation algorithms; however, discrepancies exist when one looks at the derivatives of 
the measured quantities. 



5 Numerical Results for the fractal A 

We focus now attention on the SC of fractal dimension dn = In 8/ In 3 described in section 
2 (fractal A), and provide evidence of the existence of a phase transition at a finite value 
of (3. The behavior of the SC should interpolate between the d — 1 and d = 2 cases. In 
d = 1, (3 C diverges. In fact, with a similar analysis of the partition function to the one used 
in the Appendix, it can be analytically calculated that £(/3) ~ | lntanh(/3)| _1 , which means 
that, when (3 — > oo, £(/?) ~ e 2lS . If we define (3 C {L) as the value of (3 for which £ ~ L, 
then (3 C (L) ~ InL in d = 1. This means a smooth growth; it is therefore necessary to pay 
attention in order to demonstrate that (3 C (L) does not diverge with L in the case under 
scrutiny. 

We performed different simulations for lattice sizes with L = 27, 81, 243, 729 and 2187, 
and concluded that there are several evidences of the occurrence of a phase transition in the 
Ising model on this Sierpinski carpet. Not only the average magnetization changes from 
to 1 going from the disordered to the ordered phase, but we also observed that some 
observables, shown below, present peaks increasing as powers of the size, L. The most 
stringent criterion to locate the critical coupling, j3 c is to study the crossing of the Binder 
cumulants relative to different volumes. In d — 1, where there is not a phase transition at 
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Figure 3: SDM extrapolation for the derivative of ln(|M|) in a L = 2187 lattice. Circles 
are results relative to Wolff algorithm, crosses are relative to SW. The extra point refers to 
another simulation (not a SDM extrapolation). The region where the SDM extrapolation 
has small errors is very small. 





81 


243 


729 


2187 


27 


0.6693(3) 


0.67176(5) 


0.67345(3) 


0.67455(2) 


81 




0.6730(2) 


0.67422(5) 


0.67502(1) 


243 






0.67486(7) 


0.67546(4) 


729 








0.6759(1) 



Table 1: Values of the coupling (3 where the Binder cumulants meet each other. 



a finite value of (3, the cumulants cross each other at (3 = oo. This is not the case in this 
model. The positions of the intersections of the Binder cumulants are reported in Table [1] 
and indicate that (3 C ~ 0.675. However, considering that the curves of the largest lattices 
have almost the shape of a "step-function", and that the intersections take place within 
the zone of sudden change of slope, it is difficult to obtain them with good precision. We 
have estimated the value of (3 C by other means, as we shall see in the next section. 

It has been observed |1(| that on fractal lattices, as the ones we consider, a certain num- 
ber of stages of construction of the SC is necessary in order to obtain reliable predictions 
about the critical behavior. Actually, as we shall see in the next sections, we were forced 
to discard the smallest lattice sizes, L = 27, 81, in certain fits to obtain reliable results. 



5.1 v exponent and /3 c (oo) 
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L 


# files 




0n(L) 


Max. of (ln(|M|))'(L) 


97 


~ uuuu 


u.oou 




1 ^ 1 Q M 
lo. j 


81 


' ^ (JUUU 




6975(1 ) 


28 46(6^1 


243 


~ 4000 


0.650 


0.6505(2) 


58.5(2) 


729 


~ 13900 


0.6618 


0.6622(1) 


115.4(4) 


2187 


~ 6400 


0.6682 


0.6681(2) 


218(1) 



Table 2: Maxima of (ln(|M|))'(L). 



L 


# files 


@s\ra 


Po{L) 


Max. of U'(L) 


27 


~ 6600 


0.580 


0.5744(3) 


8.40(2) 


81 


~ 5800 


0.629 


0.6252(2) 


17.56(9) 


243 


~ 4000 


0.650 


0.6493(2) 


35.8(2) 


729 


~ 13900 


0.6618 


0.6617(1) 


69.8(4) 


2187 


~ 2400 


0.668 


0.6680(3) 


132(3) 



Table 3: Maxima of U'(L). 
5.1.1 Position of the maxima 

In a FSS analysis, every observable, O, which diverges in the thermodynamic limit when 
= /3 c (oo) displays a peak for finite L, located at a different value 0o(L). When L 
increases, the difference between 0o(L) and f3 c (oo) becomes smaller according to the law 

Po{L) - &(oo) oc L" 1 /". (28) 

A three-parameter fit determines the values of exponent v and of c {po). To produce such 
a fit we have considered three different observables: the derivative of the logarithm of the 
magnetization, the derivative of the Binder cumulant (^) and the isothermal susceptibility 
([]). These three quantities display maxima at the values of /3o{L) shown in Tables 0, || 
and H, respectively. 

In these Tables (0, [| and f|) we report the values of the peaks for the different lattice 
sizes together with the statistics performed in correspondence with each L. Here and in 
the following, each file is made of 200 steps of SW, 400 steps of Metropolis and 2000 cluster 
sweeps with the Wolff algorithm. The integrated correlation time is always less than a file 
of measurements. 

The best fits are obtained discarding the L = 27 data from the three Tables. We obtain 
the results shown in Table |5|. All the fits yield % 2 /DF < 1 (where DF means the number 
of degrees of freedom in the fit); we conclude from these that: 

/3 c (oo) = 0.6751(1) v- 1 = 0.59(1), (29) 

and 

v= 1.70(1). (30) 
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L 


# files 




(3JL) 


Max. ofx(L) 


97 


^ uuuu 






97 07/"^ 


81 


~ 5800 




D 63442 f4"l 


1 79 9(2) 


243 


~ 2000 


0.651 


0.65394(7) 


1164(4) 


729 


~ 8000 


0.663 


0.66402(2) 


7759(7) 


2187 


~ 4400 


0.669 


0.66936(1) 


51883(58) 



Table 4: Maxima of x(L)- 








v- 1 


X 


0.67522(7) 


0.589(3) 


(H\M\)y 


0.6742(5) 


0.62(1) 


U' 


0.6747(7) 


0.61(2) 



Table 5: Fits to eq. ( ^8f) for different observables O, lattices 81, 243, 729 and 2187. 



We remark that we obtain a better result using the susceptibility than the other two 
observables (the derivatives of the logarithm of the absolute value of the magnetization 
and of the Binder cumulant). 



5.1.2 Fits at fixed values of the Binder cumulant 

Let us now check the values of /? c (oo) and v obtained (see Table |5|) employing eq. 
with different definitions of the apparent critical coupling (3o{L). We consider now a fixed 
value go of the Binder cumulant Ul(P), and define P(L) via the equation 

U L (P(L)) = g . (31) 



It is clear from reference |T7j that P(L) behaves like /3q(L), see equation (|28|). 

We have used for g the values 0.5, 0.6 and 0.7 (see Fig. |j). The values of the f3(L) 
given by condition ([H]) are collected in Table |6| for the lattice sizes L = 81, 243, 729, 2187, 
and the results of fitting these values to Q2"BD are shown in Table [F[ We see that these 
values of /3 c (oo) and v~ x are compatible with those given in (p9|). We remark that the 
fits reported in Table [7| have been obtained using only the lattice sizes 243, 729 and 2187. 
In all these fits the number of degrees of freedom is zero, and the reader should take the 
results of the Table |as a check of our previous estimates of (3 C and 1/u. 



5.1.3 Peaks of (ln(|M|))'(L) and U' L 

According to the finite size scaling ansatz the derivatives of (ln(|M|))(L) and Ul must 
show a peak increasing as L 1 ^. The values of the maxima in correspondence with each L 
are reported in Tables || and |[ 

In the case of U' L a two-parameter fit gives 

L > 81 : u- 1 = 0.600(4), x7 DF = 1-37, (32) 
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Figure 4: Binder cumulants for L = 81, 243, 729, 2187. 



9o 


L = 81 


L = 243 


L = 729 


L = 2187 


0.5 


0.61972(7) 


0.64677(4) 


0.66026(2) 


0.66736(3) 


0.6 


0.62551(4) 


0.64959(2) 


0.661737(6) 


0.668099(4) 


0.7 


0.63134(2) 


0.65245(2) 


0.663205(9) 


0.668868(6) 



Table 6: Values of f3(L) where the Binder cumulants take the values 0.5, 0.6 and 0.7. 



which again agrees with our previous result (|29|). 
For (ln(|M|})'(L) we obtain 

L = 729, 2187: v 



0.579(5). 



(33) 



For L > 27 or L > 81 the fits are bad (i.e. x 2 /DF is large). 

In these two cases (derivatives of U and ln(|M|)) we have not the necessary accuracy 
in order to have a stable fit taking into account the scaling corrections. 



5.2 7/V exponent: scaling of x 

The isothermal susceptibility (§) displays a peak increasing as 

By fitting the peaks with L 1 ^ we obtain a poor % 2 /DF; however, after discarding the 
data relative to L — 27, 81 we obtain a x 2 /DF < 1. The result is: 

L > 81 : j/v = 1.729(1) x7 DF = 0.6. (34) 

In order to include also the data relative to L — 27,81, we considered the scaling 
corrections to such a fit according to eq. ([24]) . Applying the procedure of the fit (eq. (|24])) 
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On 

yu 


8Joo) 


z/" 1 


0.5 
0.6 
0.7 


0.6753(2) 
0.6751(4) 
0.6752(1) 


0.585(6) 
0.589(2) 
0.584(3) 



Table 7: Fit of the values of Table | to equation (p8|) . 

to the data of Table § (i.e. we compute the effective exponent 7/z/ using all the lattice 
sizes: i.e. 7^(27, 81), 7/1/(81, 217),...) we find: 

7/z/ =1.76(2) cj = 0.27(14) x7 DF = 7, (35) 

which is unsatisfactory. However, if we discard only the L = 27 data we end up with 

L>27: 7/z/ = 1.730(1) u = 1.9(8). (36) 

In the last fit, shown in Fig. [5], we have not degrees of freedom, because we have fitted 
three data to (0), which contains three free parameters. However, the result for 7/z/ is 
compatible with (|34] ) having included the L — 81 data and the exponent u represents the 
correction to the scaling. It is clear that the final value of 7/z/ in the fit is really stable 
(almost-)independently of the u value. 




0.00000 0.00005 0.00010 0.00015 0.00030 0.00025 



Figure 5: Fit to obtain the scaling corrections in x using L = 81, 243, 729, 2187. 
We conclude that: 

cj = 1.9(8) 7/z/ = 1.730(1). (37) 
By using v = 1.70(1) we can write down the result for the 7 exponent 

7 = 2.94(2) . (38) 
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(\Mt \)(0 6751) 


(\Mr\)(0 6750) 


(|Mr|)(0 6752) 


\\ ivj -l\/ \yc\^j j 


97 










81 


W . I C/thJ I J. 1 


7Q42f 1 


7Q4QM 1 


794^^ 


243 


0.7345(1) 


0.7338(1) 


0.7351(1) 


0.7345(8) 


729 


0.6742(1) 


0.6729(1) 


0.6754(1) 


0.6742(14) 


2187 


0.61338(8) 


0.61110(8) 


0.61561(9) 


0.61336(234) 



Table 8: Absolute value of the magnetization at the critical point for the different lattice 
sizes. 

5.3 f3/v exponent: scaling of (\M\) 

Having computed the critical value of the coupling (|29]), we can extract the f3/v exponent 
from the scaling law : 

(|M L |)(/3 c (oo))ocL-^. (39) 

We have written in Table |8| the values of (|Ml|) in correspondence with the points 
0.6751, 0.6750 and 0.6752, and the final value that we take to make the fit (that is, it takes 
into account the error in our final value of (3 C ). 

Only by fitting the lattice sizes with L > 243 it is possible to obtain a "reasonable" fit 
(x 2 /DF = 2.4/1; we have only one degree of freedom), in particular 

P/ v = 0.080(1). (40) 

The quality of the data is not good for trying a fit taking into account the scaling 
corrections. 

5.4 d s and hyperscaling 

Having obtained the values of 7/V = 1.730(1) and (3/v = 0.080(1) we can consider the 
hyperscaling relation 

1- + - = d, (41) 

V V 

where d is the dimension that controls hyperscaling. In our case, we obtain d = 1.890(2), 
which is close to the Hausdorff dimension of the lattice under scrutiny. 

Finally we remark that it is clear from the standard derivation of the hyperscaling 
relations that the dimension that plays a role is the Hausdorff one. 

5.5 a/u exponent: the specific heat 

The specific heat, for positive values of a, has a peak which scales with L a l v ' . In our case, 
this peak is very flat showing a plateau that goes down as the lattice size increases. This 
corresponds to a negative value of the exponent a. Our numerical estimate is a ~ —1.15 
(which derives from vdu = 2 — a). Notice that such a value lies between the zero value of 
the two dimensional Ising model and the infinitely negative value of the one dimensional 
Ising model. 
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0c 


V 


i 


uj 


d = 1 


OO 


oo 


oo 


2 


this fractal 


0.6752(1) 


1.70(1) 


2.94 (2) 


1.9(8) 


d = 2 


0.44... 


1 


1.75 


4/3 



Table 9: Comparison between the results for the phase transition in the d = 1, d = 2 and 
d — 1.89 (fractal A) Ising models. 



L 


/^sim 


&{L) 


Max. oidM/dp 


Max. of x(L) 


64 


0.7396 


0.74548(6) 


5.797(3) 


86.57(2) 


256 


0.8060 


0.8094(3) 


9.06(2) 


839.8(5) 


1024 


0.8480 


0.8501(6) 


12.67(5) 


8484(3) 


4096 


0.8780 


0.8782(6) 


16.6(6) 


88808(101) 



Table 10: Results of the simulation of the fractal B. 



As expected the exponents for the present lattice interpolate between the corresponding 
values of the d = 1 and d = 2 shown in Table |9|. The value of the u exponent in 

d = 1 is calculated in the Appendix. 



6 Numerical results for the fractal B 

In order to provide further evidence for the hyperscaling relation we considered a second 
fractal lattice, the Sierpinski carpet B, described in section 2, whose Hausdorff dimension 
is du = In 12/ In 4. We carried out a finite size scaling analysis with lattice sizes of L = 
64, 256, 1024 and 4096. We were interested in the verification of the hyperscaling relation 



(p|). In order to get an estimate of the critical exponents j/u and f3/v from a single 
simulation for every size L, we considered in this case the observables \ and the derivative 
of the magnetization with respect to (3, as it is explained in the following paragraphs. The 



results of the simulations are shown in Table 10 



The exponent 7/1/ is obtained from the peak of the susceptibility. A fit with the three 
first lattice sizes gives the value j/u = 1.6530(2), and the value 7/1/ = 1.6755(4) with the 
three last ones. This difference means the existence of scaling corrections. We shall take 
as our best estimate the value 

7/1/ = 1.67(2). (42) 

Differentiating the magnetization with respect to (3, we obtain a quantity which scales 
as L lyl ~^ > l v . Discarding again the L = 64 data, we find: 

nn „m, .2 



0.241(3) r/DF = 3.1. (43) 

The large value of x 2 /DF again shows the existence of large scaling corrections. However, 
using only the data of the lattice sizes L = 1024 and L = 4096 we have found (1 — (3)/u — 
0.20(3) that is compatible in the error bars with (fl3|). 
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0.15 



0.20 
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L 
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Figure 6: Best fit to (3 c (po) obtained when varying u, using the values of /3 C (L) shown in 
Table 110. The error in /3 c (oo) contains the uncertainty in v. 



Now we need an estimate of the v exponent in order to get j3/u from fl^3"|). We have 
obtained a good fit to eq. (EHf) using the positions of the peaks of the derivative of M with 



respect to (3. These are shown as (3 C (L) in Table |10[ In this case we are dealing with a 
three-parameter fit, so that we need to use the four lattice sizes for consistency. Letting v 
vary freely, we obtain the best fit for 



v = 3.23(8) (3 c (oo) = 0.928(3) x7 DF = 0.95. 
This fit is shown in Fig. |]. Using this result for v and (fl3"D , we get 

P/u = 0.069(10). 



(44) 



(45) 



With the results (|42|) and fl45|), we obtain through the hyperscaling relation (|4"ID , the 
value d = 1.81(3). The error in this value is large owing to the scaling corrections, but it 
is perfectly compatible with the Hausdorff dimension of the fractal. 



7 Conclusions 

In the present paper we have considered the issue of phase transitions on discrete frac- 
tal lattices. First, by means of an exact Migdal-Kadanoff decimation transformation we 
confirmed the absence of a finite temperature critical point for the Sierpinski gasket. 

Then we focused on lattices with infinite ramification order and performed extensive 
MC simulations on deterministic Sierpinski carpets of different sizes and different Hausdorff 
dimensions, where we observed the presence of a finite temperature phase transition. The 
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FSS analysis was applied in order to extract from the data the critical exponents and their 
numerical errors. 

To the best of our knowledge the present study represents the first massive attempt to 
determine numerically the critical properties of Ising models on fractal lattices. We have 
illustrated the difficulties one encounters when dealing with these systems, which are due 
to the slow convergence of the averages to the infinite volume limit. 

Interestingly, the hyperscaling relation, which was not obtained before for such lattices, 
indicates that the role of d of the Euclidean case is taken by the Hausdorff dimension. 
Finally we should comment about the universality of the present results. As pointed out 



by Wu and Hu [|18|], on a true self similar structure the statement of universality survives 
only in a weak version. In fact, the various exponents might depend on the detailed 
structure of the fractal. It has been found, by means of approximate treatments, that the 
exponents vary even when the fractal dimension remains fixed, because they depend also 
on other geometrical characteristics such as the connectivity and the lacunarity. This is not 
too surprising since due to the scale invariance a small change at the smallest stage may 
be magnified up to the largest one. In other words, the lattice structure, which in normal 
critical phenomena does not enter with all its details here plays a much more important 
role. However, we have verified that the hyperscaling relation holds for the two fractals we 
have considered, which suggests that for lattices having the same fractal dimension, (3, v 
and 7 may vary with the lacunarity, but the sum of their ratios cannot! 

We finally remark that the model under study could provide a very good benchmark 
to study the thermodynamics of the Ising model near its lower critical dimension. 



Acknowledgments 

We wish to thank L.A. Fernandez for very useful discussions. J.M. Carmona is a Spanish 
MEC fellow. J. J. Ruiz-Lorenzo is supported by an EC HMC (ERBFMBICT950429) grant. 
The numerical work was done using the RTNN parallel machine (composed by 32 Pentium 
Pro processors) located at Zaragoza University. 



8 Appendix 



In this appendix we will compute the correction-to-scaling exponent for the one dimensional 
Ising model. 

We will use the Binder cumulant (0), which can be written as 



where 



(M n ) c 



(M 4 ) e 
2[M 2 )l 

d n In Z v 



dh n 



(46) 



(47) 



h=0 



with h the magnetic field and Zy the partition function for the volume V (in this case, 
V = L). 
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It is possible to diagonalize exactly the transfer matrix for the one dimensional model, 
and the two eigenvalues are 



A 1)2 = [cosh(/i) ± (sinh 2 (/i) + e" 4/3 ) 1/2 
Then the partition function is 



Zy — \± + \\ , 



From p7|) and fl49|) , we obtain 



(M 2 ) c = Ve 



2f3 1 - (tanh(/3)) 
1 + (tanh(/3)) 



v 



V ' 



(48) 
(49) 

(50) 



and 



(M 4 



Ve 2P 



1) (2sinh(/3)) iK - (2cosh(/3)r + 12V 2 e 4/3 (4sinh(/3) cosh(/?)) 



(2sinh(/3)) y + (2cosh(/5)) v ' 



(51) 

We put these expressions into (f5|), and consider the derivative of the Binder cumulant. 
It has a peak at finite values of (3 for each L, j3 c {L). The peak of the derivative of the 
Binder cumulant scales as L x l v plus scaling corrections. In d = 1, v = oo, so that we have 



U' L (/3 c (L)) = A + BL-" + 0(L- 



-2u>\ 



(52) 



We calculated numerically this derivative and obtained the plot shown in Fig. [7] for several 
values of L. 

Fitting to equation (|5!3) we obtain, when L — > oo, the asymptotic limit 



= 2. 



(53) 
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